Transition to the Giant Vortex State in an Harmonic Plus Quartic Trap 



H. Fu and E. Zaremba 
Department of Physics, Engineering Physics and Astronomy, 
Queen's University, Kingston, Ontario K7L 3N6, Canada. 
^ ■ (Dated: February 2, 2008) 

We consider a rapidly rotating Bose-condensed gas in an harmonic plus quartic trap. At suffi- 
£Nj ' ciently high rotation rates the condensate acquires an annular geometry with the superposition of 

a vortex lattice. With increasing rotation rate the lattice evolves into a single ring of vortices. Of 
interest is the transition from this state to the giant vortex state in which the circulation is carried 
by only a central vortex. By analyzing the Gross-Pitaevskii energy functional variationally, we have 
been able to map out the phase boundary between these two states as a function of the rotation rate 
^vq . and the various trapped gas parameters. The variational results are in good qualitative agreement 

' with those obtained by means of a direct numerical solution of the Gross-Pitaevskii equation. 
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PACS numbers: 03.75.Hh, 03.75.Lm, 67.40.Vs 



I. INTRODUCTION 



In a recent paper, Fetter et at j]| investigated the vortex structure in a rapidly rotating Bose condensate confined by 
an harmonic plus quartic trapping potential. The interest in quartic confinement stems from the fact that the rotation 
rate O is limited in harmonic traps by the trap frequency ui± , at which point the centrifugal potential destabilizes the 
system and the radius and angular momentum of the condensate diverge. The addition of the quartic potential Q 
avoids this instability and allows one to investigate rotation rates higher than loj_. In this regime, the condensate 
exhibits a more complex structure, both with regard to its density distribution and the arrangement of vortices within 

-6 : ^ mi 

The combined harmonic, centrifugal and quartic potentials give rise to a Mexican hat potential and the mean 
density acquires a local minimum on the axis of symmetry of the trap. Within a Thomas-Fermi analysis, the central 
i ^ i ■ density goes to zero at some limiting angular velocity fi^ QLIEEli beyond which the condensate takes on an annular 
structure. This behaviour is also found in numerical solutions of the Gross-Pitaevskii equation [J. Within some 

t-H , range of interaction strengths and angular velocities, these calculations reveal a structure in which a ring of vortices 
surrounds a central hole containing a multiply quantized vortex. With increasing rotation rate, a new state is favoured 
in which the annular condensate is vortex free and all the circulation is carried by a central vortex. This state with a 
multiply quantized central vortex is referred to as the giant vortex state 0, IU, El • 

In this paper we are concerned with the transition from the state with a single ring of vortices, which we also refer to 
as an annular array, to the giant vortex state. For low interaction strengths, the numerical Q, UI3 and analytical Q 
I/"") ' evidence indicates that the transition is continuous, with the radius of the ring shrinking in size until the ring is 
absorbed by the central hole. The situation at higher interaction strengths seems to be different, with the annular 
array persisting as a metastable state. In any case, one expects the giant vortex state to be preferred energetically 
above some critical angular velocity tt c . We determine the phase boundary between these two states by variationally 
minimizing the Gross-Pitaevskii (GP) energy of the annular array and comparing this with the corresponding energy 
of the giant vortex state. Our results are in essential agreement with those of Kim and Fetter Q which are based on 
the same physical model but are obtained using a different calculational approach. 

In Sec. II we present the physical model of the annular array within the context of the Gross-Pitaevskii energy 
functional. The nature of the problem suggests that a good starting point for the calculation of the GP energy can be 
achieved by replacing the annular array by a vortex sheet. Nevertheless, the discrete nature of the vortex cores leads 
to important corrections to the energy which determine the relative stability of the annular array and giant vortex 
states. We thus begin with a detailed analysis of the vortex sheet problem in Sec. Ill and then examine in Sec. IV 
the various corrections to the vortex sheet energy due to the vortex cores. Our numerical results are presented in Sec. 
V and we conclude with a discussion in Sec. VI. 

II. ENERGY 

The GP energy functional in a frame of reference rotating with angular velocity Q is given by Q 
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where is the condensate wave function, a > is the s-wave scattering length and L z = J d 3 r^!*z ■ (r x p)^f is the 
^-component of the angular momentum. The normalization of the wave function is N = J d 3 r\^\ 2 , where N is the 
total number of particles. As in previous discussions OI0,E1I3]> we consider the two-dimensional limit in which the 
confining potential is taken to be an harmonic plus quartic potential in the radial direction with no confinement in 
the z-direction: 
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Here d± = \J h/AIu>± is the harmonic oscillator length. For this two-dimensional situation the wave function W 
depends only on the radial and angular variables r and </>, respectively. 

Our purpose is to determine the equilibrium state of the system by minimizing the GP energy functional for a given 
angular velocity Q. Variations with respect to yield the GP equation which can be solved numerically to determine 
the exact equilibrium state. Solutions of this kind have been obtained Q, y, Q but the calculations are numerically 
demanding and cannot be performed in all parameter regimes. It is then useful to adopt a variational approach which 
has the added advantage of providing more physical insight into the nature of the solutions. To this end we introduce 
the amplitude-phase representation of the wave function, W = ^/crne 19 , where a is the density per unit length in the 
z-direction. The two-dimensional density n(r, <p) then has the normalization 



d 2 r n(r) = 1 . 



(3) 



Using d± as the unit of length and Hu>± as the unit of energy, the (dimensionless) energy per particle takes the form 

1 



E[n,v] = J d 2 r j^lVv^l 2 + ^nv 2 + Vn + ^grr <> r vn 



(4) 



where the condensate velocity is given by v = V0 and the interaction strength is g — A-naa. In this context the 
equilibrium state is determined by the pair of functions n and v which minimize this energy functional. Such an 
approach is feasible when, as in the present case, the physical states of interest are known. For > the condensate 
forms an annulus with some distribution of vortices P, Q . Initially the vortices are arranged on a regular lattice, but 
as the central hole becomes well-established, the lattice evolves into concentric rings of vortices around the central 
hole. Eventually, a single ring is formed and the transition from this state to the giant vortex state is the transition 
of particular interest here. The geometry and variables used to describe the annular array are illustrated in Fig. ^ 
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FIG. 1: Schematic of the annular array configuration. 



The velocity field for N r vortices on a ring of radius R is assumed to be given by a superposition of the velocity 
fields of individual vortices, 
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where Rj are the positions of the singly-quantized vortices. The second term accounts for the vorticity of a central 
vortex of strength N . In this approach, the form of the velocity field is effectively fixed, but R, N r and iVo remain 
as variational parameters. For the ring configuration of Fig. ^ the velocity and density have an angular periodicity 
of 2ir/N r , for example, 



v(r, (j) + 2n/N r ) = v(r, 0) , 
and all such quantities can be expanded in a Fourier series. We write 

v(r) = v (r) + vi(r) 
where Vq is the azimuthally averaged part of the velocity field which is given by 



vo(r) 
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This can be interpreted as the velocity field of a central vortex and a vortex sheet of radius R across which the velocity 
is discontinuous. It should be emphasized that this part of the velocity field is perfectly general and is independent 
of the assumed vortex superposition in Eq. (JSJ) . It is valid even for the exact GP solution having a net circulation of 
No in the low-density hole region and N r vortex singularities arranged on a ring of radius R. 
The remaining part vi = Vi r r + v\$(p has the Fourier expansion j^| 



N Nr (r 
V\ r (v) — sin kN r cj) I — 



kN r 
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N Nr (r \ kNr 

v uf ,(r) = sgn{r-R)—J2 coskN r<l>( — ) , (10) 

r k=l \ r > / 

where r< (r> ) is the lesser (greater) of r and R. It is this part which contains the singular behaviour of the velocity 
field at the positions of the vortices. For example, by performing the sum in Eq. (|1U|) we have 

N r /3(cos7V^-/3) 

«10 = sgn(r - R) — ——2 — — — , 11 

r 1 + p z — 2/3 cos N r <p 

where (3 = (V</r>) r . One can easily see from this expression that oc p^ 1 near each vortex, where p is the 
distance from the vortex. We also note that the vi velocity field is strongly localized near r = R since the factor (3 
decreases rapidly with increasing \r — R\ when N r is large. That is, the velocity field rapidly approaches the long-range 
part described by Vo as one moves away from the ring of vortices, and is therefore essentially azimuthal at the edges 
of the annulus. In a sense, this provides an a posteriori justification for the validity of the superposition of individual 
vortex velocity fields in Eq. Q since no geometric constraint is being imposed on the flow by the boundaries of the 
annulus. This is unlike the situation of a single vortex near the surface of a planar boundary where the velocity field 
is strongly modified from that of an isolated vortex by the presence of the boundary [Toj . 
Using Eq. J7J), the classical kinetic energy consists of the following terms: 

Ek ^\\ d2rn ^ v o( r ) + J d 2 rn(r)v {r)v lltl (r) + d 2 r n(r)«?(r) , (12) 

where 

1 [ 2 * 

n {r) = — J n(r,(p)d<f) . (13) 

We see that the separation of the velocity field into v and vi leads to a kinetic energy contribution which only involves 
vo and the angular average of the density no(r). It is thus natural to decompose the density as n(r) = no(r) + ni(r), 
where n% is the part containing all higher Fourier components. 
With these definitions, the energy functional can be written as 



E[n,v] = E vs [n , v ] + E vc [n, v] , 



(14) 
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with 



Evs[no,v ] = J d 2 r 
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and 

Ey C [n,v] = / d 2 r 



-|V\/n| 2 + -nvj + nw o wi0 - flrnv^ + -gnj 



(15) 



(16) 



The term Ey$ signifies the contribution to the energy that arises when the discrete array of vortices is replaced by a 
vortex sheet (VS) of radius R and circulation jV r . We note that the variables entering this part of the energy functional 
are independent of the angular variable </>. We identify the second term Eye as the vortex core (VC) energy. If this 
contribution to the energy is neglected, the Ey$ functional by itself is sufficient to define both no and vq. At this 
level of approximation the functional provides a Thomas-Fermi (TF) description of the vortex sheet configuration. 
We emphasize, however, that no approximations have been made in writing the energy functional as the sum of the 
two separate contributions in Eqs. i|15fl arid (|16fl . 

The second contribution Eve accounts for the actual vortex core structure and depends on the full inhomogeneous 
nature of the density and velocity field near each vortex. We note that the singular part of the velocity field arises in 
the terms containing Vi . This singularity is of course compensated by the density which must behave as n oc p 2 near 
the core of each vortex. To represent this behaviour we write the density as |ll| 

n(r) = F(r)n(r) , (17) 

where F(r) is an envelope function which accounts for the depletion of the density near each vortex with respect to 
some overall smooth background density n(r). Since n(r) is normalized to unity, n(r) in general will not be. 

The factorization of the density in Eq. I|17l) is not unique but it allows for a convenient variational representation 
of the density depletion around each vortex that occurs on a length scale £ which is small in comparison to all other 
characteristic lengths in the problem. To be specific we choose an envelope function having the form 

N r N r 

F(r) = e ~l R - Rl l 2 /« 2 - e Hr ~ Rl|2/?2 , (18) 

i=l i=l 

which gives the density depletion a gaussian profile. The parameter £ represents the vortex core radius and is treated 
as a variational parameter. It will be of the order of the local healing length and we anticipate that £ <C b, where 
b = 2nR/N r is the inter-vortex spacing. In this situation, the density close to each vortex is given to a good 
approximation by n(r) ~ (1 — exp(— p 2 / £ 2 ))n(r) , which has the required n oc p 2 behaviour. An alternative to Eq. i|18|) 
is to represent each vortex core by the piece- wise continuous function f(p) = (p/£,) 2 for p < £ and f(p) = 1 for 
P > £ 0- IS • The advantage of the gaussian core is that it allows an essentially analytic calculation of the core energy. 
However, in view of the variational nature of the calculation, one would expect the two forms to yield quantitatively 
similar results. 

Once the core structure is accounted for via the envelope function F(r), the background density n(r) is expected 
to have a weak angular variation. It can itself be treated as a variational function in the minimization of the total 
energy. However, we expect and confirm that the VC contribution to the energy is relatively small in comparison to 
Evs and as such can be treated as a perturbation. Our strategy is therefore to minimize Eys with respect to no and 
^o and to use the information provided by this minimization in the evaluation of Eye- Since the Vi velocity field 
is highly localized near r — R, the terms containing Vi in Eq. (|16[) are only sensitive to the background density in 
this region. Given its smooth angular variation, we will for simplicity choose h{r) to be a function of only the radial 
variable r. According to its definition in Eq. I|17[l . this implies that the background density is related to the vortex 
sheet density by 

n(r) = ^\, (19) 
^o(r) 

where Fo(r) is the angular average of F(r). We use this equation to estimate n near r = R. In this way all 
contributions to Eye become functions of £ and to complete the calculation, the VC energy is minimized with respect 
to this parameter. 
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III. THE VORTEX SHEET APPROXIMATION 



In this section we calculate the energy within the vortex sheet approximation (VSA) obtained by minimizing Eq. I|15|) 
with respect to the density n and the circulation parameters of the azimuthally averaged velocity given by Eq. JHJ . 
Here we denote the circulation parameters by v\ — N n and v 2 = N n + N r . It is reasonable to treat these parameters 
as continuous if the circulation is large. To impose the normalization constraint on the density we minimize the free 
energy Fyg — Eyg — jj, j d 2 rn where \i is the chemical potential. Variation of Fyg with respect to n gives 

, , / Mi ~ Ux{x), r < R , , 

9Mr) = { , 2 ~ U 2 \x{ r>R, ( 2 °) 

where 

fii = /i + VLVi , (21) 

and 

U i (x) = ^(^+x + \x^J (22) 

with x — r 2 . The inner and outer radii, Ri = ^fxl, of the annulus are defined by 

Vi = Ui(xi), » = 1,2. (23) 

By its definition, uq is the azimuthally averaged density and is therefore a continuous function of the radial variable 
r. The requirement that uq be continuous across the vortex sheet leads to the relation 

*>-i(3 + 3). (-) 

which indicates that the average of the velocity on either side of the vortex sheet is equal to the velocity for rigid 
body rotation at the radius r = R. Eq. I|24|) can also be interpreted in terms of the circulation provided by a uniform 
vortex lattice. The number of vortices contained within a circle of radius R that is needed to give the velocity R£l at 
r = R is h*L = i? 2 fi. Thus, Eq. Q24J) is equivalent to the relation vl = [v\ + v^) /2. 

The density given by Eq. 1)20(1 also depends on the chemical potential [i. For a given v\ and vi, this parameter is 
fixed by the normalization of the density in Eq. This gives the additional relation 

dx{fi! - Ui(x)) +7T / dx (fi 2 - U 2 (x)) , (25) 

J Xo 



where xq = R 2 . 

Using Eq. J5UJ), the free energy can be expressed as 



F VS = -y / dx (/i! - U^x)) 2 -y dx{^ 2 - U 2 (x)) 2 . (26) 

The equilibrium state corresponds to the minimization of Fys with respect to the two parameters v\ and v 2 . Noting 
that Xq, X\ and X2 are implicitly functions of v\ and 1/2, this variation yields the equations 

f x ° r° dx 

fl dx( fH -U 1 (x)) = v 1 — - U^x)) , (27) 

Jx± Jx± % 

n / dx ( M2 - U 2 {x)) = u 2 — 0*2 - U 2 (x)) . (28) 



The six nonlinear equations, Eqs. (|23I25|) . (|27|l and (|28|) . are sufficient to determine the six unknown parameters v\, 
V2, R, Ri, R2 and [i. 

The procedure described corresponds to a free variation of the circulation parameters v\ and i>2- Alternatively, the 
minimization can be carried out with imposed constraints, such as fixing the radius R of the sheet or the number of 
vortices N r in the ring. In the first case, the pair of equations (|27|l and 1)28(1 is replaced by the single equation (R 
fixed) 



r° r 2 r 2 dx 

n dx(in-u x (x))-n dx(fi 2 - u 2 {xj) + 2x n / —(^-u 2 {x)) 

J X\ J Xq J 



X 
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dx 



(Mi - Ui(x)) + 



dx 



(M2 - U 2 (x)) 



(29) 



while in the second (N r fixed) we have 



dx (m - Ui(x)) +Q dx(n2- U 2 {x)) - N r 



dx 



(M2 - U 2 (x)) 



dx 



(Hi - Ui(x)) + 



dx 



(M2 - U 2 (x)) 



(30) 



In the limit N r — > 0, (z^i = ^2) the latter equation reduces to the minimization equation for the giant vortex state: 



n 



dx 

dx (/ii - Ui{x)) = ui I — - Ui(x)) 



(31) 



The parameter xo of course has no significance in this limit. We later make use of some of these constrained 
minimizations. 

Although the above equations can be solved numerically (to be described later), it is useful to perform a perturbative 
analysis in order to gain some insight into the form of the solutions. For Q > fi^,, the condensate is confined to an 
annulus whose radius increases with Q. At the same time, the width of the annulus decreases. We therefore expect 
the parameter w — x 2 — x\ to be small in comparison to xq. In this situation, we can find useful relations between 
the various parameters by expanding the integrals in Eq. I|27[l in a Taylor series, 



f x ° 1 1 

/ dx f(x) = f(x 1 )A 1 + -f'(x 1 )A 2 1 + -f"( Xl )Al + 
J X1 2 6 



(32) 



where Ai = xq — X\. In our case, f{x\) — and the expansion starts with terms of order Af . The expansion of the 
integrals in Eq. 1271) can thus be viewed as providing a power series expansion for v\. 



V\ — v x + v 1 Ai + ■ 



(33) 



We insert this expansion for v\ wherever it appears in Eq. (|27|l . including the /ii and U\(x) terms, and generate a 
power series in Ai for each side of the equation. Equating the coefficients of like powers of Ai and retaining terms to 



order A*, we find 



V\ = flxi 



2Ai a x /A 



where 



Analyzing Eq. 128|) in a similar way, we obtain 

v 2 = flx 2 



1 



3 si 18 V x\ 



1 - 2n 2 + Xx! 
1 - n 2 + 2\xi 



2A 2 a 2 /A 



(34) 



(35) 



(36) 



3 x 2 18 V x 2 

where A 2 = x 2 — x and a 2 is obtained from the expression for ai by replacing X\ by x 2 . These relations show that 



1 



Xq 



2Q 



(vi + v 2 ) = - (xi+ x 2 ) 



1 { a l A 2 , "2 A 2 



12 \ Xl 



x 2 



(37) 



that is, R 2 = (R 2 + R 2 )/2 with corrections of order Af . These equations also imply that Ai ~ w/2 to lowest order. 

To obtain an explicit expression for x\ + x 2 we take the difference of the two equations in Eq. (|23|l . thereby 
eliminating [x. Using the above expansions in powers of u>, we then find x\ + x 2 ~ (O 2 — 1)/A, which implies 



x 



n 2 - 1 
2A ' 



(38) 
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with corrections of order (w/xo) 2 . This leading order result for xq is the same as the expression for (xi+x 2 )/2 obtained 
for a vortex lattice with a hole 1]. The reason for this can be seen by inserting the lowest order result vi = VlXi 
into Eq. J23). This gives fi = \ (xi(l — Q 2 ) + Xx 2 ) which is recognized as the equation giving the boundaries of 
the annulus in the vortex lattice state. This correspondence indicates that the annular array in the largc-f2 limit is 
effectively undergoing rigid body rotation. 

The result for xq in Eq. 1[38[) shows that the denominator of ol\ is 1 — fl 2 + 2\x\ ~ —Aw. Thus the terms formally 
of order Af in Eqs. (|34|1 and l|36(l in fact give rise to corrections to Vi that are of order w. From Eqs. (|34|) and 
we thus find 



N r 



V 2 - V\ 



1 



i n 2 - 1/3 
I n 2 - i 



(39) 



Another useful relation follows from the expansion of the normalization condition, Eq. i|25|) . To order w 3 we have 



7T 

9 = o 



' dU 2 




dUi 




2 T 


' d 2 u 2 




d 2 U x 




dx 


X2 


dx 


Xl_ 


w 

48 


dx 2 


H 

X2 


dx 2 





w 



(40) 



The first term on the right hand side appears to be of order w 2 , but if the potential derivatives in this term are 
evaluated to lowest order in w we obtain X(x 2 — xi), which in fact is of order w. Thus, to include all contributions of 
order w 3 , the potential derivatives in the first term must be expanded to first order in w and the second term must 
also be retained. Doing so one finds 



I2g 



1/3 



n(X + n 2 /2x ) 
For large O we have w ~ (Gg/nX) 1 ^ 3 = w^. Since w 

d = R 2 — Ri = 



I2g 



1/3 



X Wqq 

2~n 



nX(i + n 2 /(n 2 - 1)), 

R\ , the physical width of the annulus is given by 



1 

3^2 



(41) 



(42) 



that is, d oc il 1 for large fl whereas the radius R of the vortex sheet is proportional to fl. We also note that 
the strength of the central vortex is Nq — flxo — N r /2 ~ f2 3 /2A for large S7, while that of the vortex sheet is 

N r ~ 5^00/12. 

To obtain accurate results for all O we solve Eqs. I|23I25|I . (|27|l and (|28|l numerically. We have found that this can be 
done in a straightforward iterative manner. We start by using the large-il expressions v\ = S7 3 /2 A, v 2 = ^i+5f2w oc /12, 
xq = (il 2 — 1)/2A and x 2 = xq + Woo/2. We then estimate an approximate chemical potential as \i! = U 2 (x 2 ) — Qv 2 , 

define /ii = /i' + flvi and determine x\ from U\(x\) = fi\. We now calculate the integrals Fi(a, b) — J dx (/ii — Ui{x)) 

and Gi{a,b) = J^dx(/ii — Ui(x))/x in order to determine new values of the circulations from Eqs. I|27|l and l|28|l : 
v\ = flFi(xi,xo)/Gi(xi, xq) and v 2 = D,F 2 (xq,x 2 )/G 2 (xo,x 2 ). At the same time we use Eq. (|25H to determine the 
effective coupling strength g' = n(Fi(xi, xq) + -F^ieoj 2^2))- The chemical potential which gives the desired value of g 
is fx = fx' + A/i and the change in chemical potential is estimated as Afi = (g — g')/n(x 2 — x\). With these updated 
values of v\ , v 2 and /i we recalculate xq — [v\ + is 2 )/2tt and Xi from U{xi) — fii, and repeat the remaining steps 
described. Convergence to an accuracy of one part in 10 6 typically required 10-20 iterations. 
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477 


52.2 


7.42 


8.44 


7.93 


1.02 


0.954 


0.174 



TABLE I: Physical parameters within the vortex sheet approximation for A = 1/2 and g = 1000. The values in parentheses 
are obtained from the solution of the GP equation The last column gives the vortex core radii as determined in Sec. IV. 
All lengths are in units of d±. 

In Fig. |21 we show the equilibrium density n$(r) within the VSA for g = 1000 for a few values of fl. The density 
exhibits a downward cusp at r — R which corresponds to the average depletion with respect to a smooth background 
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FIG. 2: The solid curves show the density as a function of r within the vortex sheet approximation for three values of Q; 
A = 1/2 and g = 1000. The dashed curves are the densities for the giant vortex state. All quantities are in units of 



that the vortex cores give rise to. To make this more apparent, we have also plotted the density for the giant vortex 
state which only contains the central vortex. This shows the effect of transferring circulation from the core of the 
giant vortex to the vortex sheet near the middle of the annulus. The significant depletion of the density that occurs 
at the vortex sheet radius is compensated by an increase in the width of the annulus. The fact that the depletion 
exhibits a cusp is of course an artifact of treating the vortex sheet within a Thomas-Fermi-like approximation. In 
reality, the cusp will be smoothed out on a length scale £ characterizing the size of the vortex cores. Nevertheless, 
the overall qualitative behaviour is expected to be a good representation of the azimuthally averaged density in an 
annulus containing a ring of vortices. We note further that there is a relatively slow recovery of the density from the 
region of the sheet to the outer regions of the annulus. This is associated with the long-range behaviour of the vortex 
sheet velocity field. For a single vortex in a uniform gas of density n^jy, the size of the core is set by the healing length 
£o = V \/87ra«3D, but the density actually approaches its asymptotic value quite slowly: n(r) — ri3£>(l — £o/?" 2 + ■•■)• 
A similar behaviour is occurring in the present context, but it is partially masked by the finite width of the annulus. 
Another feature of interest is the nearly constant maximum value of the density as a function of Q. We can define 
an average density n in the annulus by 2nRdn — 1, where d is the width of the annulus. Since R is approximately 
proportional to £1 while d is inversely proportional to f2, the near constancy of n, and hence the maximum, follows. 
The same argument also applies to the giant vortex state. 

The various parameters that emerge for the VS state are collected in Table I. The asymptotic dependences R ~ 
\J (SI 2 — 1)/2A, Nq — f2 3 /2A and N r ~ ^Iw^lYl are found to work quite well down to Q = 4. Also given in 
parentheses in the table are the values of the parameters as determined by a numerical solution of the Gross-Pitaevskii 
(GP) equation We see that the VS approximation underestimates the circulation of the GP solution by about 
10% while the discrepancy in the radius is only a few percent. The discrepancy in the width d of the annulus is much 
larger, with the VS width approximately 30% smaller than the GP result. It should be noted, however, that there is 
some arbitrariness in the way that the GP widths are extracted from the numerical data since the GP densities fall 
off smoothly to zero [l2T |. In fact, despite its inherent uncertainty, visual inspection of the published figures yields 
widths which are in much closer agreement with our VS results. We thus believe that the difference in the widths is 
actually much smaller than indicated in Table I. Finally, the difference between our inter-vortex spacings b and the 
GP results simply reflects the different values of the circulations N r in the two calculations. Interestingly, we see the 
same slight increase of b with increasing Q as seen previously . We conclude that the VS approximation provides a 
reasonably faithful average representation of the state containing a ring of vortices. 

In Fig.|21we show Eys and Egv as a function of for g = 1000 and A = 1/2. To leading order in f2, both of these 
energies behave as — f2 4 /8A. The difference between them can be seen to be small and is displayed in more detail 
in Fig.0] We see that Eqv — Eys — 4fio;_L over the range of f2 shown, that is, the vortex sheet always has a lower 
energy than the giant vortex. This is to be expected since the insertion of a ring of vortices brings the velocity field 
closer to that of rigid body rotation. Likewise, the insertion of another vortex sheet would be expected to lower the 
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FIG. 3: The solid curve gives the vortex sheet energy Evs and the dashed curve the giant vortex energy Egv (in units of 
hu>±), as a function of the angular velocity fi (in units of co±), for g — 1000 and A = 1/2. 
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FIG. 4: The solid curve is the difference, Bgv — £Vs, between the giant vortex and vortex sheet energies, as a function of the 
angular velocity f2, for g — 1000 and A = 1/2. The dashed curve is the vortex core energy Eye- The point of intersection of 
these two curves defines the critical angular velocity f2 c beyond which the giant vortex state is more stable. 



energy further. As we shall see, this is indeed the case. However, the relative stability of these vortex sheet states is 
overstated since the energy associated with the vortex cores has not been included. The vortex core energy will at 
some point make states with a large number of vortices less stable. In the following section we calculate the vortex 
core energy to correct the energy as determined in the VS approximation. 
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IV. CORE ENERGY 



We begin by considering the first term in Eq. JTSJl which represents the quantum mechanical kinetic energy of the 
condensate. The main contribution to this integral comes from the region of the vortex cores which are defined by 
the envelope function F(r). Using Eq. (|17fl and treating h as a smooth function, the first term in Eq. l|16f) gives 



E 



(i) 
vc 



d 2 r 



n|VF| 2 
8F 



N r 



d 2 r 



h\\7F\ 2 
8F 



(43) 



where we have used the periodicity of F to reduce the integral to the area A defined by < r < oo and —n/N r < cf> < 
ir/N r . We assume that a vortex is positioned at r = R, = within this cell. If the core size £ is small on the scale 
of the inter-vortex spacing, as turns out to be the case, the envelope function can be approximated within the cell as 



F(r) ~ 1 - e 



_ P - P 2 /e 



(44) 



where p is the distance from the vortex centre. Assuming also that £ is small on the scale of the variations of n, we 
find 



where d = ir(ir 2 /6 - l)/2 = 1.01306- 
solely through h(R). 

The next contribution to Eye is 



Using Eqs. © and (|rU|l . we have 



E% ~ C x N r n{R) , 



(45) 



It should be noted that the dependence of this contribution on £ arises 



p(2) _ 



d 2 rnv 



N r n(R) / d 2 r{l~e- p2 ' e )vl 



(46) 



ft 2 



l-(3 2 



1 + 2 P k cos kN r (j) 



k=l 



l + f3 2 - 2f3cosN r 



(47) 



The latter form is useful in order to see that v\ oc p 2 near each vortex. This singularity is cancelled by the envelope 
factor in Eq. (|46|l . However, in order to evaluate Ey^, the first line in Eq. (|47|) proves to be more useful despite the 
appearance of the factor (1 — /? 2 ) -1 which is singular at r = R. As we shall see, this singularity is removable. 
To proceed we simplify the envelope factor using p = |r - R| = ^Jr 2 + R 2 -2ri?cos0 (see Fig.QJ: 



!_ e -p7r = 1 _ e -(r 2 +R 2 )/e e 2rRcos4>/e 
~ 1 _ e -{r-Rf/f e -rR<f/e 

= (1 - e -(-«) 2 /« 2 ) + e ~(r~R) 2 /e (i _ e -^ 2 /? 2 ) . 

In going to the second line we have made use of the fact that the angular range in is small when N r is large while the 
grouping of terms in the last line is introduced to facilitate the cancellation of singularities that arise. Accordingly, 
we consider the integral 



' A 

2ttN,. 



d 2 r{\ - e~ 
dr 



-Rf/e 



i 



-(r-Rf/e 



h 



il 2 



p 2 



(48) 



It is apparent that the singularity at r — R is cancelled by the zero in the envelope factor. We show in Appendix A 
that this integral can be evaluated analytically to a good approximation, with the result 
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where M 2 = |(iV r - ±)(N r - §) and a = \e*l 2 = 0.6673..., with 7 being Euler's constant 13]. The logarithm is 
associated with the (r — singularity of the integrand in Eq. 148|l which is cut off at the distance £. The appearance 
of N r in the logarithm indicates that £/& is the physically relevant parameter, with b the inter-vortex spacing. This 

is also true of the dominant power law terms in the expansion which is given to higher order in Appendix A. 

(2) 

The second contribution to Ey C involves the integral 



d r e 



(1-e- 



Nz 



l-p 2 



w/N r 
-ir/Nr 



1-e 



-rR<j> 2 /C 



l + 2^/3 fc coskNrC, 



k=l 



(50) 



In this case all Fourier components of v\ contribute. We require the discrete Fourier transform of a gaussian which is 
given by 



r/iVr 



2irg m = N r e~ ri ^ ^ cosmN r <t>i 



-ir/N r 



-m 2 /4/t 



K[erf(z m )] 



(51) 



where k = rR/£ N 2 , z m = tt^/k + im/2y/~K and erf (2) is the error function [13j. The inverse Fourier transform is 



E 



Qrn 6 



imcj) 



(52) 



from which follows the useful identity 



E 9m = 1 



(53) 



m— — 00 



With these results the angular integral in Eq. (|50(l can be expressed as 



1-e" 



■n/N r 



) 



1 + 2 P m cosmNr 



2tt - J^erf^^) - 2 J /3 m e - m2/4K 5R[erf (z m )] 
V K V K TO=1 

— oc 

-^(l-/3 m )e-™ 2 / 4K 3?[erf(z m )]. 

m=l 



Thus we find 



2N r 



(54) 



We now see that the (1 — /3)" 1 singularity is cancelled by the numerator, leaving the factor J2T=o / ( 1 + /3) . Both 
the gaussian and (3 2 factors are highly localized around r = R, which allows the sum to be evaluated at r = R to a 
good approximation. We thus find 



where 



JO r 
1 — 00 

C 2 = ^ me- m2 / 4K 5R[erf(^ m )] . 



(55) 



(56) 



The integrals in Eq. (|55(l on the ranges (0 < r < R) and (R < r < 00) are examples of the integrals defined in 

Appendix B. The constant C 2 can be expressed conveniently as the integral 



C 2 = 2k 



)cot — e 



(57) 
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which shows that C2 behaves asymptotically for large k as C2 — 2 v / 7tk(1 — + • • •) = 2y/irR/£ i N r + ■ ■■. Since 

the integral in Eq. 1)55(1 is proportional to £ as £ — > 0, I2 itself has a finite limiting value in this limit, unlike the 

(2) 

logarithmically divergent behaviour of I\ . Our final result for E vc is 



E^ c = ^N r h(R)(h + h). (58) 



This quantity is plotted as a function of £ in Fig. 03 
The next contribution to the core energy is 



E VC = / d 2 rnv v 1(i> 

~ -JV r n(i2) / d 2 re-^ 2 ^e- rR ^^v v^. (59) 



This is the cross term that arises from squaring Vo + vi. For r > R, vo and Vx<p have the same sign, indicating that the 
velocity is actually larger than Vq close to the sheet. On the other hand, has the opposite sign for r < R, implying 
that there is a cancellation between vq and v\$ on this side of the sheet. The overall negative sign in Eq. (|59[1 is due 
to the fact that only the regions within the vortex cores contribute to the integral. In these regions the density is 
depleted with respect to the average no and Ey C compensates for the contribution these regions make to the VSA 
kinetic energy. 



The angular integral in Eq. (|59)l is the same as encountered previously in the calculation of Ey\, . We have 



E% = -N r n(R) / dr e~ {r ~ R)2 ''' e u (r)sgn(r ~ R)J- V f3 m e - m2/4K di[eri(z m )} . (60) 

Jo V « m=1 

We retain only the rapidly varying j3 m factor from the sum in doing the radial integral and set r — R in the remaining 
slowly varying factors. We then obtain 

— 00 

- J2 [(No + N r )J> Nr+1 (R/0 - NoJ^^R/O] e-" l2 / 4K 5R[erf(z m )] , (61) 
hi 

m=l 

where the integrals are defined in Appendix B. Within the sum there is a strong cancellation between the 

two Nq terms, leaving mainly the term proportional to N r . This is most evident in the £ — > limit for which the 
integrals take a limiting value of y/ir£ t /2R. The remaining sum is evaluated using Eq. I|53|) . giving Ey^ c ~ 
—N r n(R)ir 3 / 2 (N r t;/2R) to lowest order in £. This shows that Ey^ depends linearly on £ for £ — > 0. 

The contribution Ey C — — J cPrQrnvic/) is essentially of the same form as Ey C and can be handled in a similar 
way. We find 

— 00 

- E [ J ^N r -iW0 - J^N r+ xW0] e- m2 / 4 ^[erf(z m )] . (62) 
hi 

m=l 

We recall that Q.R 2 = Nq + N r /2, with Nq typically being much larger than N r . However, the cancellation between 
JmN -1 an d JmN +i within the sum diminishes the size of this contribution in comparison to EyQ, as can be seen in 

Fig. in 

Finally, we have the interaction term 

El = l 9 J d2rn2 

= l -N r g [ d 2 r[l-h(r)] 2 h 2 (r) 



2 

= \g J d 2 m 2 + \g J d 2 r[(h 2 ) - (h ) 2 ]n 2 (r) , (63) 

where h(r) = exp (— (r — i?) 2 /£ 2 ) exp (— rR<fr 2 /t; 2 ) and the subscript in the second term indicates the zeroth Fourier 
component as defined in Eq. (|51J) . The first term in Eq. 1)63(1 contributes to Eys, while the second term represents 
the interaction contribution to Eye- We thus find 



= $9 d 2 r[(h 2 ) -(h ) 2 ]n 2 (r) 
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FIG. 5: The total vortex core energy in units of hu)±_ as a function of the core radius £ in units of d±, for Q = 6, g = 1000 
and A = 1/2. The numerically labelled curves indicate the various components Ey^, i = 1, • • • , 5. The dashed curve gives the 

(2) (3) 

combined result E vc + E V q, as this represents the total change in classical kinetic energy with respect to the VS state. The 
minimum in the total core energy gives the optimal core radius. 

In Fig. we show the various contributions to the core energy as a function of £, together with the total Eye — 
Si=i^vc' for the case g = 1000 and fl — 6. The behaviour for other values of these parameters is very similar. 

(3) 

All the contributions except for E vc are seen to be positive. As explained previously, the latter is the correction 

(2) (3) 

to the kinetic energy that arises from the interference between vq and x>i$. Since E vc and E vc both represent 
corrections to the VS kinetic energy, we have also plotted the sum of these two terms. The combination is seen to 
decrease monotonically with increasing £. The other contributions Ey C , Ey C and Ey C are positive and increase 
with £. The total core energy Eye exhibits a minimum at £ m i n ~ 0.179 which is the equilibrium core radius for 
the case being considered. Values of £ m i n for other are given in Table I. The core radius is seen to be rather 
insensitive to the rotation rate. These values are of the order of, but somewhat larger than, the local bulk healing 
length £ ~ l/^8Traan(R) ~ 0.12. 

The dependence on £ of the quantum kinetic energy term Ey^ is entirely due to the n(R) = n Q (R) / F (R) factor, 
which also appears in all the other core energy contributions. The angular average of the envelope function is given 
by Fo(R) = 1 — er{(w^/K)/2^/^^K and its dependence on £ appears through k = (R/N r £) 2 . For k > 1 the error function 
is approximately unity and F (i?) _1 ~ 1 + v / 7r(^/6) + 7r(£/6) 2 + • • • where b = 2-KR/N r is the inter- vortex spacing. At 
the minimum, £ m in/fr — 0.2 and the Fo(R)^ 1 factor makes h(R) about 50% larger than no(R). This is perhaps a slight 
overestimate of the background density since the gaussian envelope function gives a core density profile that is more 
compact than what one would expect it to actually be. As explained earlier, the core density recovers its asymptotic 
value rather slowly due to the slow decay of the azimuthally averaged velocity. This latter effect is accounted for in 
an average way within the VSA but is not accounted for using the gaussian core profile. In Fig. we compare no(r) 
and n(r); the enhancement of n(r) above no(r) is evident. However the detailed structure in n(r) should not be taken 
seriously. The cusp at r = R is of course a residual artifact of the cusp in no(r); it would be eliminated if the smooth 
core profile were accounted for in no(r). Apart from this, the figure suggests that fi(r) would appear smoother if the 
gaussian where augmented by wings which decayed more slowly so as to be more consistent with the long-range part 
of the core density profile contained in no(r). To give an impression of what an improved n(r) might look like, we 
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FIG. 6: The solid curve shows the background density, h(r) = no(r)/Fo(r), as a function of r for fi = 6, g = 1000 and 
A = 1/2. The dashed curve is the density in the vortex sheet approximation, no(r). The chain curve is a possible refined 
smooth background as discussed in the text. 

have plotted the function nfit(r) = A(r — i?i) (i?2 — r) with A chosen to reproduce n(R). With this choice of A one can 
see that hfit(r) goes smoothly to no(r) at the edges of the annulus and appears to be a reasonable candidate for the 
true smooth background. The ratio no (r) /fist (f) can be thought of as a refined Fq(t) which in turn implies a refined 
core density profile. It should be emphasized, however, that the detailed shape of n(r) is not particularly relevant 
since only n(i?) is used in estimating the vortex core energy. In view of the variational nature of the calculation, 
refinements in the core density profile would not be expected to lead to significant changes in the value of the core 
energy. 

V. RESULTS 

We now combine the core energy with the VS energy calculated in the Sec. III. To begin, we do this perturbatively, 
that is, we assume that the core energy is a small correction to the VS energy. This is certainly correct for large but 
will become less accurate as fl is reduced. Nevertheless, for the time being we use the parameters R, N r and no(-R) 
which arise from the minimization of Eys to evaluate Eye over the whole range of Q. The procedure provides an 
upper bound to the total energy Eys+Eyc of the annular array and therefore underestimates its stability relative to 
the giant vortex state. To determine the transition point we plot Eqy—Evs and Eye versus fi. This is done in Fig.0] 
for the case of g = 1000. As stated earlier, Eqv—Evs is almost constant as a function of f2, while Eye increases 
approximately linearly. This latter behaviour is due to the fact that all the vortex core energies are proportional to 
N r which itself is approximately a linear function of il. It should be noted that the vortex core contributions also 
depend on the parameter N r £/ R, but this parameter is approximately constant since R is also proportional to O. 

The intersection in Fig. 21 occurs at the angular velocity fl c ~ 7.1. Repeating the calculations as a function of 
g yields the phase boundary between the vortex array and giant vortex states. These results are presented as the 
dash-dot curve in Fig. [7| and are very similar to those obtained by Kim and Fetter p| as indicated by the open 
squares. These authors do not invoke the VS approximation but rather represent the density as in Eq. 1|17[) . choosing 
the envelope function to be composed of linear cores and the background density h(r) to have an analytic form that 
was found to work well for the giant vortex. They then fix the number of vortices in the ring, N r , and minimized 
the GP energy with respect to the remaining parameters. They find a local minimum in this energy as a function 
of N r which corresponds to the annular array and then compare this energy with that of the giant vortex (N r = 0) . 
Remarkably, the minimum values of N r determined by this approach agree very well with our VS results in Table [fl 
even though the latter excludes our core correction. However it should be emphasized that the VS approximation 
already includes the vortex cores in an average sense, which is why our VC energy is typically a relatively small 
correction. For g = 1000 and A = 1/2 they find f2 c ~ 6.7, slightly smaller than our perturbative result of 7.1. The 
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FIG. 7: The dash-dot curve shows the phase boundary between the annular array (AA) and the giant vortex (GV) state as 
determined by treating the core energy perturbatively. The solid points joined by solid lines denote the position of the phase 
boundaryas determined by a global minimization of the total energy. The open triangles are the values determined by the GP 
solution [l| (the g = 1000 point is a lower bound for fl c ) and the open squares are the results of Kim and Fetter ||. The filled 
square with error bars gives the approximate bounds on Q, c as determined in Ref. Q for g — 125. The dashed curve denotes 
the onset of a hole (VLH) in the vortex lattice (VL). 

fact that our critical angular velocity is slightly higher would suggest that our variational treatment of the annular 
array is somewhat better, but given the differences in the approaches, the two results should be considered as being 
essentially in agreement. Their results in fact approach ours for smaller values of g. 
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FIG. 8: Total energy contours (Evs + Eve) in the vicinity of the global minimum (GM). The point labelled VS is the location 
of the minimum in the vortex sheet approximation (Evs)- The parameters in the calculation are £7 = 5, g = 1000 and A = 1/2. 

A more accurate calculation would involve minimizing the total vortex array energy (Eys + Eye) with respect 
to the two parameters v\ and v-i , or alternatively, N r = — v\ and R = (y\ + i^)/2f2. In Fig. [S] we present a 
contour plot showing the behaviour of the total energy in the vicinity of the global minimum. It can be seen that 
N r and R are almost orthogonal variables. The point labelled VS is the position of the minimum of the vortex sheet 
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energy Ey$ itself; constant energy contours of Eys would appear similar to those shown but would be centred on 
VS. As one moves away from VS along a line at constant R, Eyg of course increases but the total energy actually 
decreases as a result of the decrease in Eye- The fact that Eye decreases with N r is somewhat surprising since Eye 
is explicitly proportional to N r . However, Eye also depends on the parameter £/&, and the inter-vortex spacing is 
decreasing with increasing N r . This dependence is evidently dominating the variation of Eye with N r . 

To locate the global minimum we have found it convenient to begin by using the fixed- N r minimization scheme 
in Eq. I|30[l which determines the minimum value of Eys for a given value of N r . R is found to vary only slightly 
as N r is varied and the minimum of the total energy along this line in the N r -R plane can readily be determined. 
With N r fixed at the value where this minimum occurs, we subsequently evaluate the total energy as a function of 
R, adjusting the chemical potential to ensure that Eq. I|25ll is satisfied. The global minimum can then be approached 
quickly with successive variations of N r and R. In this way we find a global minimum at N r — 37.5 and R = 4.83 for 

= 5. These values are perhaps coincidentally close to the GP results of 37 and 4.83, respectively; the agreement 
with the GP values (in brackets) is slightly poorer for f2 = 6 where we find N r — 47.0 (44) and R — 5.87 (5.80) and for 

= 7, N r = 57.5 (51) and R = 6.88 (6.85). Nevertheless, the global minimization in all cases moves N r in the right 
direction from the VS values in Table I, and yields values of R that are in excellent agreement with the GP results. 
It is therefore clear that our variational approach is providing a good description of the annular array properties at 
this value of g. 

By repeating these calculations as a function of Q and comparing the energy at the global minimum with Egv-, 
we obtain an improved value of Q c . These points are plotted in Fig. and show that the phase boundary shifts to 
slightly larger angular velocities; the fractional change in f2 c is 3-6% and increases with decreasing g. The phase 
boundary for the global minimum ends at g = 200; below this interaction strength the total energy did not exhibit a 
local minimum corresponding to the annular array. 
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FIG. 9: The density in the vortex sheet approximation for (N r = 10, No =4), solid curve, and (N r — 12, JVo = 1), dashed 
curve. The results are obtained for fl = 2.5, A = 1/2 and g — 125. 

We believe this apparent failure of the calculations at low values of g is associated with the breakdown of the 
Thomas- Fermi approximation used for the vortex sheet. To confirm this we have analyzed in more detail the case of 
g = 125 at ft = 2.5 which was studied previously by Kasamatsu et al. Here, we calculated the total energy as a 
function of integral values of Nq and N r . With N r fixed at 10, we find a minimum energy at Nq = 4 which matches 
the (N r = 10, No = 4) equilibrium state found in Ref. The VS density for this state shown in Fig. [5] corresponds 
quite well to the density given in Fig. 1 of Ref. Q with regard to the boundaries of the annulus, the radius of the 
vortex array and the relative density of the inner and outer portions of the annulus. However, the total energy we 
find is higher than that of the giant vortex state (—1.35 hu>± vs. —1.41 hoj±); according to our calculations, £1 = 2.5 
lies to the right of the phase boundary, consistent with the extrapolation of the dot-dash curve in Fig. but not 
consistent with the stability of the (10, 4) state as determined by the GP equation. As stated earlier, we attribute 
this discrepancy in f2 c to the breakdown of the TF approximation for the VS. In fact we find states with even lower 
total energy when the central circulation is minimized. For N r = 12 we find a minimum of E = —1.5Huj± at Nq = 1 
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and the corresponding density is shown by the dashed curve in Fig. [5J It is clear from this figure why the energy 
of this state is lower in our approach. Although the VS energy is higher than for the (10, 4) state, the vortex core 
energy Eye of the (12, 1) state is reduced to an even greater extent because of the low density no(R) at the radius 
of the vortex sheet. (Recall that the VC energy is proportional to this quantity.) However we expect that there will 
be significant beyond-TF corrections to the VS energy for this state due to the highly inhomogeneous nature of the 
density. We believe this is the reason why the (10, 4) state, as opposed to the (12, 1) state, is the lowest energy 
state as found in the GP calculations Q|. We conclude that beyond-TF corrections must be playing a role in our 
determination of the phase boundary for small values of g. 
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FIG. 10: The density in the vortex sheet approximation for one (solid), two (dash) and three (dash-dot) vortex sheets. The 
results are obtained for fl — 5, A = 1/2 and g — 1000. 

We finally present some results for configurations involving multiple rings of vortices. These results are obtained 
using a straightforward generalization of the vortex sheet approximation discussed in Sec. 1 1 1 1 1 for a single ring. In 
Fig. 1101 we show the density no (r) obtained by minimizing the VS energy functional for one, two and three vortex 
sheets for the case of ft = 5. The energies for these cases are —127.0, —128.2 and —128.7 hus±, respectively, indicating 
that the vortex sheet approximation prefers multiple ring configurations. In fact, one can consider the limit of a 
continuous distribution of vortex sheets with a circulation density v{r) which is a continuous function of r. The 
distribution which minimizes the free energy is v(r) — Q/tt, that is, a constant density throughout the annular region 
of the condensate. This gives rise to the rigid body velocity field ^rb (t) — rfl as used in the uniform vortex lattice 
calculations Of course the vortex core energies must be taken into account in order to determine the relative 
stability of the multiple vortex sheet states. When this energy is included we expect to find a sequence of phase 
boundaries between the n and (n+ 1) sheet states, and in particular, a phase boundary between the single and double 
ring configurations which lies to the left of the phase boundary shown in Fig. \7\ Our calculation of the vortex core 
energy can be extended to treat these multiple ring configurations but we will not pursue this extension here. 

VI. CONCLUSIONS 

In summary, we have investigated the properties of an annular BEC in an harmonic plus quartic trap in the regime 
of high angular velocities Q > luj_. Of particular interest is the transition from the state containing an annular array 
of vortices to the giant vortex state in which the circulation is carried by a single central vortex. The phase boundary 
defining the transition between these two states was determined by a variational analysis of the Gross-Pitaevskii energy 
functional. Our approach identified two contributions to the total energy. One, based on an azimuthal average of the 
density and velocity, defined what we refer to as the vortex sheet approximation. This part accounted for most of the 
energy but important corrections to the energy coming from the vortex cores had to be included in order to properly 
describe the transition to the giant vortex state. The vortex core energy was also treated variationally by assuming 
the core to have a gaussian density profile. With this choice of the core profile we were able to provide analytic 
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expressions for the various contributions to this part of the energy. The phase boundary determined by treating the 
core energy perturbatively differed only slightly from the phase boundary determined by a global minimization of the 
total energy. However, the global minimization did provide equilibrium parameters which were in better agreement 
with those obtained from the solution of the GP equation [l[ . Our results for the phase boundary are expected to be 
less reliable in the limit of weak interactions (g < 250) where corrections to our Thomas-Fermi-like approximation to 
the vortex sheet appear to be important. Our approach can also be used to deal with multiple ring configurations 
and possibly with other arrangements of vortices in rotating trapped gases. 

As a final comment, we emphasize that we have only addressed the static equilibrium properties of the vortex 
structures considered. The dynamic stability of these states is also of interest and was addressed in the paper by Kim 
and Fetter Their analysis leads to the conclusion that the vortices in the annular array are indeed dynamically 
stable. Thus the annular array is a well defined equilibrium state that is separated by a meaningful phase boundary 
from the giant vortex state. 



APPENDIX A 



We present here details of the calculation leading to the result given in Eq. I|49(l . The contribution to the integral 
in Eq. (|48|l from the range R < r < oo can be written as 



I> = 2irN r 



dx 



fl_ e -WC) 2 (*- if) 1 

V J X™r - 1 



(Al) 



Integrating by parts we have 



It 



2nR 2 



2 r°° 



dx(x - i) e -(fl/f) 2 C*-i) 2 l n 



~2N r 



x 2N r _ l 



(A2) 



The x T factor in the logarithm gives the integral (with x = 1 + y) 



A = 



2irN r £ 
R 



dye- R2y2/e y\n{l + y) 
2 2R 4 \RJ 2\R 



The other logarithmic term is expanded as 



]n(x™ r - 1) ~ H2N r y) + M t y + M 2 y 2 + 



(A3) 



and gives the contribution 



B = — 



2nR 2 



Similarly, the range < r < R gives the integral 



dye- R2y2/ ^y \\n{2N r y) + M lV + M 2 y 2 + 



It 



2, T .V,. /^(l-e-^ i' J " 



U2 /■! 



r 2N r 



1 - X 2N r 



dx{l-x)e-( R ltf^-^\n(\-x 2N r) 



2vri? 

'~? Jo 

o n>2 poo 

-—pr~ \ dy e- R y /« y [\n(2N r y) - M lV + M 2 y 2 + 
s Jo 



(A4) 



where the assumption that R/£ is large allows the upper limit to be extended to infinity. Adding I± to !{" = A + B 
we finally obtain 



h = 2tt 



N r £J 2 R 



M 2 + 



(A5) 



where a = ±e 7/2 with 7 equal to Euler's constant, M 2 = ±(N r - |)(JV r - §) and M 4 = ~j^{N r - ^)(8N? + AN 2 
218A,. + 251). This expansion provides an accurate representation of 1% over the range of £/R of interest. 
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APPENDIX B 

In Sec. IIVI we had integrals of the form 

/ d ^-*?/e {R/r) ^) = / e -^/e {1+y) -* dy ^j> m) . (Bl) 

JR r JO 

For large R/£ we have to a good approximation 

poo 

JkWO = / e- R2 y 2 /t 2 e- kU ^ dy 
Jo 

poo 

j o 




-erfc(^/a-y-) 



where a_ = (R/0 2 - fc / 2 an d y_ = k/2a-. 
The other integral that appears is 

[ R ± e -^ 2 ^ (r/R)^ = f ^-^-1^ (1 - yf dy = J< (R/0 . (B2) 
Jo r Jo 

For large R/£ this becomes 

J<(R/t) = [ 1 e-W?e-'>W-v)dy 
Jo 



/ e -R 2 v 2 /e e -k{vWm dy 

Jo 




where oq_ = {R/£,) 2 + an d 2/+ = k/2a + . 
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